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ANALYSIS OF UNSTEADY THERMAL BOUNDARY LAYERS 

by Roy W. Miller 
Lewis Research Center 

SUMMARY 

The velocity and temperature distribution in an unsteady thermal boundary layer is 
analyzed by an approximate integral method. The case of a flat plate in a free stream 
having small harmonic velocity oscillations about a steady mean value is treated in de- 
tail. The resulting velocity profiles agree with the available experimental data. These 
profiles are used to predict the thermal convection when the plate and the stream are at 
constant but different temperature levels. The flow oscillations cause harmonic tem- 
perature oscillations plus a small steady feedback effect which reduces the mean heat 
transfer to the plate. The amount of decrease is represented by a single coefficient, 
which is presented graphically. 


INTRODUCTION 

Unsteady thermal boundary layers occur in many applications of unsteady flight, 
rotating blade motions, and nozzle flows. Previous analytical investigations of unsteady 
thermal boundary layers have been asymptotic in nature and thus apply to restricted 
ranges of interest (refs. 1 to 8). For flow with an oscillating free stream, solutions are 
asympotically valid for very low or very high frequency fluctuations. Results are gen- 
erally limited to surface heat transfer, and detailed temperature distributions are not 
readily available. In this work, a general method of analysis is presented which will 
allow the calculation of approximate temperature profiles for a class of laminar two- 
dimensional flows with a wide range of applicability . 

An integral method of analysis for unsteady laminar boundary -layer flow was re- 
cently investigated (refs. 9 and 10). That method basically consists of the following 
procedure: 

(1) Obtain asymptotic solutions for limiting cases of the flow under consideration. 

(2) Assume velocity profiles with sufficient generality to include the asymptotic 
solutions. 



(3) Apply integral relations (velocity -weighted averages of the momentum equation) 
and compatibility conditions (normal derivatives of the momentum equation evaluated at 
the boundary) to determine the unknown velocity profile form parameters. 

The present method of analysis for thermal boundary layers follows a similar procedure 
by utilizing the thermal energy equation as follows: 

(1) Choose temperature profiles with sufficient generality to include known asymp- 
totic solutions. 

(2) Determine profile form parameters to satisfy temperature -weighted averages 
of the thermal energy equation. 

The method is applied to the case of a flat plate in a free stream having small har- 
monic velocity oscillations about a steady mean value. Velocity profiles obtained by the 
method of references 9 and 10 are used as input for the thermal analysis. Those results 
were obtained by using integral relations and compatability conditions. In an attempt to 
improve the accuracy of previous results, the velocity profiles have been recomputed, 
using integral relations only, with results of the new computations presented herein. 


INTEGRAL ANALYSIS OF UNSTEADY THERMAL BOUNDARY LAYERS 

The integral method of analysis for unsteady thermal boundary layers is formulated 
here for incompressible flow over a flat plate with time -dependent free-stream velocity. 
In the analysis, it is assumed that the plate and the free stream are at constant but dif- 
ferent temperatures. By establishing the general features of the method, other effects, 
such as prescribed heat flux, compressibility, or spatially varying pressure gradient, 
can later be included in the development. 


U = Uoo(t) 

T * I,*, - constant 


V(X Y,t) 

Y 

t ' TfX.Y.t) 1 — UtX.Y.t) 



Figure 1. - Flat-plate boundary-layer configuration and nomenclature. 
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Boundary-Layer Equations 


We consider the laminar two-dimensional incompressible flow and heat transfer for 
the flat-plate configuration shown in figure 1 . The constant -temperature free stream 
has velocity which depends only on time and flows parallel to the isothermal plate. If 
viscous energy dissipation is neglected, the boundary -layer equations are 


3U + 9V = 0 

3X 3Y 

2£+u^+v2H=i/ _ 1 9? 

at 3X 3Y 0Y 2 p 3X 


(1) 

( 2 ) 



(3) 

(4) 


(Symbols are defined in appendix A.) The boundary conditions are as follows: 


U(X,0,t) = 0 (5) 

V(X,0,t) = 0 (6) 

T(X,0,t)=T w (7) 

U(X,Y,t) — U^t) as Y — (8) 

T^Y^-T, as Y-co (9) 


Since equation (3) implies constant pressure in the normal direction, the imposed pres- 
sure gradient is determined by evaluating equation (2) at the free stream, namely, 



Equation (10), of course, applies for free-stream variations which depend on time alone. 
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Integral Relations 


In references 9 and 10, a system of integral relations is derived for the unsteady 
boundary -layer flow. The momentum equation (2) is multiplied by U k (k is a non- 
negative integer) and integrated in the Y-direction from the plate surface to the free 
stream. That is, a set of velocity -weighted averages of the momentum equation are 
formed. The resulting flow integral relations are as follows (see appendix B): 



For k = 0, equation (11) gives the familiar momentum integral equation. The second 
equation with k = 1 is a mechanical energy balance for the boundary layer . No such 
interpretation can be made for k > 1 . 

A set of integral relations for the temperature can be obtained in a manner similar 
to the preceding by considering the thermal energy equation. First, introduce the di- 
mensionless temperature ratio as follows: 


e = 


T - T 

w 



( 12 ) 


With this definition, the values of dimensionless temperature at the wall and free stream 
are zero and 1, respectively. We next multiply equation (4) by 6 and integrate in the 
Y-direction from the plate surface to the free stream. The result is 
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Upon imposing the flow continuity requirement and integrating by parts, we obtain the fol- 
lowing set of temperature -weighted averages of the thermal energy equations (see ap- 
pendix C) : 



The first of these equations (k = 0) is the so-called heat-flux equation (refs. 11 and 12). 
The remaining equations (k > 0) have no known physical significance. 


Method of Solution 

The integral method described in references 9 and 10 utilizes the integral relations, 
equation (11), together with a set of compatibility conditions to determine approximate 
velocity profiles. Compatibility conditions are normal (Y) derivatives of the momentum 
equation (2), with the result evaluated at the plate surface. The U-velocity component is 
taken to vary exponentially in the Y-direction, and the X-dependence is included through 
velocity profile form parameters. The exponential form obviates the definition of a dis- 
crete boundary -layer thickness. The variable form parameters are required to satisfy 
the integral relations and the compatibility conditions. The integral relations demand 
that the profiles approximate the velocity values in some average sense across the 
boundary layer and that the compatibility conditions impose restrictions on the slope of 
the profiles at the wall. 

The results in references 9 and 10 for oscillating flow over a flat plate show excel- 
lent agreement with experimental data (ref. 13) for the first-order unsteady effect. The 
second -order mean-flow profile approximation, at small values of the oscillation fre- 
quency, agrees only qualitatively with the exact asymptotic solution. Since the loss of 
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accuracy may have been caused by the nature of the compatibility conditions, the analy- 
sis presented herein employs the integral relations alone to determine the profile param- 
eters. 

The thermal boundary layer is computed in a manner similar to the flow. The as- 
sumed temperature profile parameters are determined in such a way that the thermal 
integral relations, equation (14), are satisfied. 


FLOW OVER A FLAT PLATE WITH AN OSCILLATING FREE STREAM 

In this section we conduct the detailed integral analysis of the flat-plate boundary 
layer for a free stream which performs small harmonic flow oscillations about a constant 
mean value. The flow is determined independently of heat transfer; in fact, the approxi- 
mate solutions for the velocity boundary layer are later used in the thermal analysis. 


Integral Relations 

For the case of a free stream with small harmonic oscillations, of frequency w and 
amplitude cUq, about a steady mean value Uq, the free-stream velocity in equation (8) 
is expressed in complex form as follows: 

Uco(t) = U Q (1 + ee iwt ) (15) 

The mathematical form of equation (15) suggests expansion of the local velocity in 

powers of the small amplitude parameter e. Of course, the leading term (order e®) is 

the steady -flow boundary -layer velocity with free-stream velocity Uq. The e* term 

corresponds to fluctuations, at frequency w, which vary in amplitude and phase through- 

o 

out the layer. Due to the quadratic nature of the boundary -layer equations, the e 
terms consist of oscillations at twice the fundamental frequency plus a steady feedback of 
the first-order oscillations on the mean flow (cf. ref. 10). Thus, we expand the 
U-velocity component, in dimensionless form, as follows: 

u = u Q (x,y) eu 1 (x,y)e iwt + e 2 [u 2 (x,y) + u 2 (x,y,t)] + ^(e 3 ) (16) 


where 
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( 17 ) 


In equation (16), is the complex velocity component parallel to the plate with asso- 
ciated fluctuation magnitude and phase angle. The velocity component u 2 is the contri- 
bution of the oscillations to the time -averaged flow. 

We substitute equation (15) and the expansion, equation (16), into the flow integral 
relations, equation (11) in normalized form, and obtain the following set of integral 
relations: 

Order e® : 
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Order e 1 : 



Order 


(time- independent part): 


_d_ 

dx 


l [ h **-* 


k + 2 ; 
k + 1 


k+1 



r 



(20) 



(k - 1 ) (k - 2) 
4 
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where denotes the complex conjugate of u^. 

We proceed with the integral analysis by assuming approximate forms for the veloc- 
ity component (uq, u^, and u^) and then determining the form factors to satisfy the pre- 
ceding integral relations . 


Basic Mean Flow 

It was mentioned earlier that the zeroth-order term in the velocity expansion (16) is 
the steady-flow boundary -layer velocity with constant free-stream velocity Uq. For the 
flat plate in a uniform stream, this basic mean flow is the familiar Blasius velocity pro- 
file Ug. This is, of course, a similarity solution in terms of the similarity variable, 
namely , 

V = -2— (21) 

Thus, the velocity component along the plate is 

u B = Fg(rj) (22) 

where Fg is the well-known Blasius stream function. 

Note, from equations (19) and (20), that the zeroth-order velocity appears in the 
calculation of higher order terms. The use of tabulated values of the Blasius profile is 
inconvenient since the subsequent integral relations are algebraic in nature. Thus, we 
will compute algebraic approximations for the Blasius profile by the integral method of 
analysis. 

Following references 9 and 10, consider four -parameter profile functions of the 
form 


u Q = 1 - e" A?? (l + Bt? + Crj 2 + D?] 3 ) 

These functions have the following desired behavior: 

Uq = 0 at y = 0 

> 

u 0 -1 as y-« > 


(23) 


(24) 
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Four conditions are required to determine the four constants (A, B, C, and D). They are 
the four integral equations (18) with k = 0, 1, 2, and 3. 

First, note the property 


d_ 

d? 


J* GO -OO 

/ n(7j)dy = / Sl(rj)dri 

•'0 Jo 


where 



Then, the four integral relations, equation (18) reduce to the following: 


(25) 


(26) 



- l)dr; = -2 



(27) 



Substituting the assumed profile forms, equation (23), into these relations and carrying 
out the indicated operations yields four algebraic equations for A , B , C , and D . 

It is convenient to represent the profile approximation, and its derivative, in the 
following manner: 
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(31) 


u 


0 





(32) 


where 


f. 

i 



£ Nm "™" 1 

m=l 


(33) 


In equation (33) we have introduced the following arrays: 




h im ~ 


1 

1 

A - B 


s i = 


0 

B 

AB - 2C 


0 

C 

AC - 3D 


0 

D 

ADI 


> 


(34) 




By using these representations, integral terms of the type which appear in equations (27) 
to (30) can be given function definitions as follows: 


FM) 



i 1+ 3 


-2 - (s I +s J ) 




We recall the following integration formula : 



(35) 


(36) 
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where n is a nonnegative integer and Cq > 0. Then, the function definitions reduce 
to the following : 


F(I,J) 


hlihj^i + j-2)l] 

(=1 - 


Sum on i and j from 1 to 4 


( 37 ) 


The F notation is readily extended for additional subscripts, and the following simpli- 


fied form for equations (27) to (30) is obtained: 

[F(2, 1) - F(2, 2)] - 2(A - B) = 0 (38) 

[2F(2, 1) - 3F(2, 2) + F(2, 2, 2)] - 4[F(3, 3)] = 0 (39) 

[3F(2, 1) - 6F(2, 2) + 4F(2,2,2) - F(2,2,2,2)] - 12[F(3, 3) - F(3,3,2)] = 0 (40) 

[4F(2, 1) - 10F(2, 2) + 10F(2, 2, 2) - 5F(2, 2, 2, 2) + F(2, 2, 2, 2, 2)] 

- 24[F(3, 3) - 2F(3, 3, 2) + F(3, 3, 2, 2)] - 0 (41) 


The form factors were computed to satisfy equations (38) to (41) by applying the 
computational scheme outlined in appendix D. The resulting values are as follows: 


A" 


"1.9059 

B 


1 . 5709 

C 


0.8159 

_D_ 


1.1655 


Values of the corresponding velocity profile are presented in table I, together with the 
Blasius profile values and the reference 9 results based on three integral relations plus 
one compatibility condition. The present values are closer to the Blasius data, except 
near the plate surface. 


Velocity Oscillations 

We compute the first-order oscillating component of velocity by following the proce- 
dure of references 9 and 10: that is, obtain asymptotic solutions for the flow under 
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TABLE I. - COMPARISON OF FOUR -PARAMETER PROFILE 


APPROXIMATIONS WITH THE BLASIUS PROFILE 


Simi- 
larity 
variable, 
v = y/C 

Blasius 

value 

Three integral relations 
(ref. 9 result) 

Four integral relations 

Value 

Percent 

error 

Value 

Percent 

error 

0.20 

0.06641 

0.06579 

-0.94 

0.07368 

10.95 

.40 

.13277 

.12888 

-2.93 

.14456 

8.88 

.60 

.19894 

.19097 

-4.00 

. 20710 

4.10 

: 80 

. 26471 

.25414 

-3.99 

.26518 

.18 

1.00 

.32979 

.31923 

-3.20 

.32312 

-2.02 

1.20 

0.39378 

0.38584 

-2.02 

0.38311 

-2.71 

1.40 

.45627 

.45273 

-.78 

.44526 

-2.41 

1.60 

.51676 

.51833 

.30 

. 50833 

-1.63 

1.80 

. 57477 

.58112 

1.10 

.57054 

-.74 

2.00 

.62977 

.63985 

1.60 

.63014 

.06 

2.20 

0.68132 

0.69362 

1.81 

0.68566 

0.64 

2.40 

.72899 

.74194 

1.78 

.73613 

.98 

2.60 

.77246 

.78462 

1.57 

.78099 

1.10 

2.80 

.81152 

.82176 

1.26 

.82010 

1.06 

3.00 

.84605 

.85364 

.90 

.85364 

.90 

3.20 

0.87609 

0.88068 

0.52 

0.88195 

0.67 

3.40 

.90177 

.90337 

.18 

.90555 

.42 

3.60 

.92333 

.92223 

-.12 

.92498 

.18 

3.80 

.94112 

.93776 

-.36 

.94082 

-.03 

4.00 

.95552 

.95045 

-.53 

.95360 

-.20 

4.20 

0.96696 

0.96075 

-0.64 

0.96383 

-0.32 

4.40 

.97587 

.96905 

-.70 

.97195 

-.40 

4.60 

.98269 

.97570 

-.71 

.97836 

-.44 

4.80 

.98779 

.98099 

-.69 

.98338 

-.45 

5.00 

.99155 

.98519 

-.64 

.98729 

-.43 

5.20 

0.99425 

0.98850 

-0. 58 

0.99031 

-0.40 

5.40 

.99616 

.99110 

-. 51 

.99265 

-.36 

5.60 

.99748 

.99313 

-.44 

.99444 

-.30 

5.80 

.99838 

.99472 

-.37 

.99581 

-.26 

6.00 

.99898 

.99595 

-.30 

.99685 

-.21 

6.20 

0.99937 

0.99690 

-0.25 

0.99764 

-0.17 

6.40 

.99961 

.99763 

-.20 

.99824 

-.14 

6.60 

.99977 

.99820 

-.16 

.99868 

-.11 

6.80 

.99987 

.99863 

-.12 

.99902 

-.08 

7.00 

.99992 

.99896 

-.10 

.99927 

-.06 




consideration; assume profiles with sufficient generality to include the asymptotic solu- 
tions; and determine the unknown profile form parameters to satisfy the flow integral 
relations . 

Qua si -steady flow . - The asymptotic solution of interest for the plate in an oscil- 
lating stream is the so-called "quasi-steady" flow. For the limiting case of very low- 
frequency fluctuations, the boundary -layer response is quasi-steady at the prevailing 
value of the free-stream velocity. The solution then follows the Blasius form based on 
the instantaneous free-stream velocity, that is, 



We expand equation (43) in a Taylor series for small free-stream perturbations and 
identify the first-order quasi-steady velocity as 

U 1Q = + - t?F£(t7) (44) 

2 

An approximation to the quasi-steady flow velocity is obtained by the integral method 
by first assuming the four -parameter profile form as follows: 

u^ = 1 - e -a77 (l +bt] + crp + d^**) (45) 

The constants a, b, c, and d are obtained by demanding that u^ satisfy the four in- 
tegral relations, equations (19) with k = 0, 1, 2, and 3. 

Since the qua si -steady velocity depends on the similarity variable, equations (19) 
will reduce to algebraic relations in a manner similar to the Blasius flow approxima- 
tion. We extend the previous matrix definitions to include the quasi-steady approxima- 
tion and its derivative, that is, 


u lq 




(46) 


du lq 

dr7 



(47) 


and 
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( 48 ) 


o' 


1 

0 

0 

0 " 

A 


1 

B 

C 

D 

A 

h. = 
1 m 

A - B 

AB - 2C 

AC - 3D 

AD 

a 


1 

b 

c 

d 

a 


a - b 

ab - 2c 

ac - 3d 

ad 


Then, the first-order flow integral relations can be reduced to the following: 

[F(4) + F(2) - 2F(4, 2)] - 2(a - b) = 0 (49) 

[2F(4) + 4F(2) - 6F(4, 2) - 3F(2, 2) + 3F(4, 2, 2)] - 8[F(5, 3)] = 0 (50) 

[3F(4) + 9F(2) - 12F(4, 2) - 12F(2, 2) + 12F(4, 2, 2) + 4F(2, 2, 2) - 4F(4, 2, 2, 2)] 

- 12[2F(5, 3) + F(3, 3) - 2F(5, 3, 2) - F(4, 3, 3)] = 0 (51) 

[4F(4) + 16F(2) - 20F(4, 2) - 30F(2, 2) + 30F(4, 2, 2) + 20F(2, 2, 2) - 20F(4, 2, 2, 2) 

- 5F(2, 2, 2, 2) + 5F(4, 2, 2, 2, 2)] - 48[F(5, 3) + F(3, 3) - 2F(5, 3, 2) 

- F(4, 3, 3) - F(3, 3, 2) + F(5, 3, 2, 2) + F(4, 3, 2, 2) ] = 0 (52) 

The profile form parameters were computed from equations (49) to (52) by the 
method of appendix D with the following result: 

a 
b 
c 
d 


1.40244 

0.91137 

0.68585 

■0.53889 


(53) 


The associated profile is shown in figure 2. Also shown in figure 2 are the exact quasi- 
steady solution and the reference 9 results obtained by using three integral relations plus 
one compatibility condition. Some improvement was realized by applying the fourth inte- 
gral relation. 

Stokes flow solution . - Another asymptotic solution, which will later be of interest, 


15 




Figure 2 . - Four-parameter profile approximations for first-order quasi- 
steady flow. 


is the Stokes flow solution for large values of £ = yoXX/Ug. For this limiting case, 
viscous and unsteady effects dominate and the flow depends only on y. The well-known 
Stokes flow result is the following (cf. ref. 14): 


u 


IS 


= 1 - 



(54) 


Profile factors . - The general first-order velocity profile approximation, which ap- 
plies for all frequency values, is assumed to be of the following form: 

Uj = 1 - e -Q,y (l + /3y + yy 2 + 6y 3 ) (55) 


where the complex profile form factors (a, /3, y, and 6) are functions of the variable 
? = ^coX/Uq. Equation (55) shows the desired behavior, which is zero value at the plate 
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surface and approaching unity at the free stream. Of course, for lower frequencies the 
independent variables combine as the similarity variable and equation (55) yields the 
quasi-steady approximation. Also the Stokes flow, equation (54) is a specialized form of 
equation (55) . 

We determine the form factors in equation (55) to satisfy the integral conditions (19) 
The integration over y is carried out as done previously, by introducing appropriate 
function definitions. However, since the general dependence of u^ on x is not known, 
the integral relations will result in a set of first-order ordinary differential equations for 
a, /3, y, and 6. The resulting set of nonlinear equations is of the following standard 
form: 


p. d® + p. d£ + p. Qv. + p. d5 - R_ 
11 d£ i2 d£ i3 d? i4 d? 1 


i = 1, 2, 3, 4 


(56) 


The coefficients in equation (56) depend on a, /3, y , 6, and £. Refer to appendix E for 
the detailed coefficients. 

Since equations (56) are linear in the form parameter slopes, we can solve for these 
values by applying Cramer's rule and obtain the following results: 


da 

d? 




7r 1 (a,j3,y,6,0 


d£ 

d? 


ff 2 (a,/3,y,6,C) 


dy 

d? 


> 


v 3 (a, 0 ,r ,6,?) 


(57) 


d5 

dC 


7r 4 (a;,/3,y,6,0 

J 


These complex equations are equivalent to eight real first-order differential equa- 
tions; and thus, eight initial values are required to obtain the solution. The eight values 
are the real and imaginary parts of a, /3, y, and 6 evaluated at the low-frequency 
limit (i.e. , £ — 0). The real parts of the form parameters are associated with the quasi- 
steady approximation (eq. (45)), whence 
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as £ — 0 


(58) 


a r 


An 

— i 

^r 

/V 

bA 

?r 


cA 2 

_ 5 r 


_dA 3 _ 


In fact, the quasi-steady behavior can be identified as the leading real term of an expan- 
sion for small £, namely, 


“q ~|( 1 +a l? + a 2^ + ' 

< 3 q~!( 1+b l« + ¥ 2+ ' 
1 'q~^( 1+c l ?+c 2 ?2 + 
5 q ~ ^ + <¥ + d 2? 2 + 

r 


) + ia?^l +a x ? +a 2 ? 2 + . . .) 
) + ib?(l + bjS +b 2 ? 2 + . . .) 

+ ic ^1 + Cjt; + 'c 2 £ 2 + • • 

.) + i|(i + d 1? + d/ + . . .) 


(59) 


Now, the initial behavior of the imaginary parts of the form parameters can be obtained 
by substituting the expansions (59) into equations (56) and retaining the dominant ima- 
ginary terms. Upon solving the resulting equations by the method of appendix D, we ob- 
tain the following: 


~ a i 


"0.32372 f 

(*i 

/V 

-0.51773 ? 

A 


0.07723 

6. 

L 


-0.51783 A 


(60) 


The set of equations (57) was solved numerically by using the Runge-Kutta formulas 
listed in appendix F. The initial values were obtained by evaluating equations (58) 
and (60) at £ = 0.01 - a starting value for which the quasi-steady solution was repro- 
duced. Results for the real and imaginary parts of the form parameters are shown 
graphically in figures 3 and 4, respectively. The computation was terminated when sub- 
sequent calculation showed that the Stokes flow had been approximately achieved. 
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Figure 3. - Real parts of first-order velocity profile form factors. Figure 4. - Imaginary parts of first-order velocity profile form factors. 




Profile results . - After the profile-form factors have been determined, the first- 
order velocity profiles are readily computed by using equation (55). We denote the mag- 
nitude of fluctuations by |u^ | , and a positive phase angle represents a phase lead with 
respect to the free-stream variation. 

The development of the first-order velocity for small values of £ = ^ cvX /Uq is 
shown in figures 5 and 6. Results are plotted against the similarity variable 77 to show 
the departure from quasi-steady flow. Recall that the quasi-steady fluctuations are "in 
phase" throughout the flow field. The unsteady acceleration force applies simultaneously 
throughout the flow field. The slower moving fluid near the wall responds to this force 
faster than the main stream, thus causing a phase advance near the wall. 

Profiles of the first-order velocity for large values of £ = f/wX/Uo are presented 

in figures 7 and 8. The variation as a function of y is of interest here; and, results 

o 

show that the Stokes flow limit is approximately approached at £ =3. 

Computed first-order velocity results are compared with the data of Hill and 
Stenning (ref. 13) in figures 9 and 10. Again, for small values of coX/Uq, results are 

plotted against rj ; and, for large values of coX/Uq, the normal coordinate y/Vi is 
used. In all cases the agreement between present and experimental results is very good. 

The frictional drag is related to the normal derivative of the velocity, with the re- 
sult evaluated at the wall, that is, 




e iwt ♦ 


(61) 


The first-order term again has components in phase and out of phase with the free- 
stream oscillations. Also, since the first-order term is oscillatory, there is no net 
first-order contribution to the drag in the time- average sense. 

We define a first-order wall shear stress coefficient, which is normalized by the 
magnitude of the Blasius shear, as follows: 



(62) 
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Figure 12. - Phase angle of first-order wall shear stress coefficient. 

Values of the shear coefficient components and the associated phase angle $ T are shown 
in figures 11 and 12. The present results span the asymptotic limits given by Illingworth 
in reference 2. 

Effect of Oscillations on the Time-Average Flow 

We complete the flow analysis by computing the steady second-order velocity pro- 
file, that is, the feedback of the first-order oscillations on the steady mean flow. Again 
we consider a limiting-case solution and then apply the integral method of analysis. 

Quasi -steady flow . - The asymptotic solution for quasi-steady flow (low-frequency 
limit) is obtained by continuing the expansion of equation (43). The resulting second- 
order quasi-steady velocity is as follows: 



This solution is indicated by the solid-line curve in figure 13. 

We consider approximate quasi-steady profile forms as follows: 

u 2q = e" p?7 (qr 7 + r rj 2 + stj 3 ) (64) 
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Figure 13. - Second-order quasi-steady velocity profile approximation. 


These forms have zero values at the wall and free stream, and they admit an inter- 
mediate zero for certain values of the coefficients. 

The constant profile factors in equation (64) are determined in a straightforward 
manner by the integral method. We first expand the array definitions (48) as follows: 
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Then when the F-function notation is introduced, the four integral relations with 
k = 0, 1, 2, and 3, equations (20), are reduced to 

[F(4, 4) - F(4, 1)] + 2[F(6, 1) - 2F(6, 2)] + 4q - 0 (66) 

[3F(4, 4) - 4F(4, 1) + 6F(4, 2) - 3F(4, 4, 2) - 2F(2, 1)] + 2[F(6, 1) - 6F(6, 2) 

+ 3F(6, 2, 2) + 8F(7, 3) + 2F(5, 5)] = 0 (67) 

3[2F(2, 2) - 3F(2, 1) - 3F(4, 1) + 8F(4, 2) + 2F(4, 4) - 4F(4, 4, 2) - 4F(4, 2, 2) 

+ 2F(4, 4, 2, 2)] + 12[4F(7, 3) + F(5, 5) - 4F(7, 3, 2) - F(5, 5, 2) 

+ 2F(5, 3) - 2F(5, 4, 3) + 2F(6, 3,3)] + 2[12F(6, 2, 2) 

- 12F(6, 2) - 3F(6, 1) - 4F(6, 2, 2, 2)] = 0 (68) 

5[3F(2, 2) - F(2, 2, 2) + 6F(4, 2) - 6F(4, 2, 2) + 2F(4, 2, 2, 2) + F(4, 4) - 3F(4, 4, 2) 

+ 3F(4, 4, 2, 2) - F(4, 4, 2, 2, 2)] - 8F(4, 1) + 4F(6, 1) - 12F(2, 1) + 12[4F(7, 3) 

+ F(5, 5)] + 5[F(6, 2, 2, 2, 2) - 4F(6, 2, 2, 2) + 6F(6, 2, 2) - 4F(6, 2)] 

+ 12[4F(7, 3, 2, 2) + F(5, 5, 2, 2) + F(3, 3) - 2F(4, 3, 3) + F(4, 4, 3, 3) 

- 8F(7, 3,2)- 2F(5, 5, 2) ] + 48[F(5, 3) - F(5, 4, 3) - F(5, 3, 2) 

+ F(5, 4, 3, 2) + F(6, 3, 3) - F(6, 3, 3, 2)] = 0 (69) 

where the property expressed by equation (25) has been observed. 

The solution of equations (66) to (69) is obtained as outlined in appendix D, with the 
result as follows: 


p 


1. 06634 ~ 

q 


0.063332 

r 


0. 161027 

_s_ 


-0. 056376_ 


(70) 
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Values of the corresponding quasi-steady profile approximation are shown in figure 13. 
The result of reference 9 is also shown for a five -parameter profile with coefficients 
determined by applying three integral and two compatibility conditions. The use of only 
integral relations gives a better approximation, even though the number of profile pa- 
rameters is one less. 

Profile factors . - We generalize equation (64), to apply for all frequency values, as 
follows: 


u 2 = e~ <py (t^y + py 2 + Xy 3 ) 


(71) 


The profile-form factors (<p, \p, p , and X) are, of course, functions of the variable 

? = J/coX/Uq. We identify the limiting behavior of the functions as the quasi-steady re- 
sult, namely, 


’< 


'pA' 


rsj 

qA 

"q 


rA 2 



sA 3 


as £ -0 


(72) 


The computation of the form factors follows the first-order flow analysis: that is, 
substitute the profile forms into equation (20) with k = 0, 1, 2, and 3; introduce the 
appropriate function notation and carry out the integration over y; and rewrite the re- 
sulting differential equations in the standard form. The result is 


P il~ +P i2~ +P i3^ 


+ P i4 


dX 

d? 


= R. 


i = 1. 2, 3, 4 


(73) 


where coefficients and R^ are listed in appendix E . 

The first-order differential equations (73), with the initial conditions (72), were 
solved numerically; the results are shown in figure 14. The associated velocity profile 
development is shown in figures 15 and 16 with t] and y, respectively, as independent 
variables. The profiles show spreading of the effected zone into the stream. 

The steady second-order contribution to the wall shear stress is related to the 
y -derivative of the u 2 profile, with the result evaluated at the wall. We define a shear 
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coefficient as the ratio of wall shear stress to the exact quasi-steady value, namely, 



( 74 ) 


Values of the shear coefficient are compared with previous analytical results in figure 17. 



u 0 


Figure 17. - Second-order wall shear stress coefficient. 


THERMAL BOUNDARY LAYER ON A FLAT PLATE WITH 
AN OSCILLATING FREE STREAM 

The thermal boundary layer is computed for the flat plate in a free stream with 
small harmonic velocity oscillations about a steady mean value. The free stream and 
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plate are at constant but different temperature levels. As a result of the velocity oscil- 
lations, the boundary -layer temperature performs small oscillations about the basic 
mean temperature distribution. Also, there is a small steady feedback of the oscilla- 
tions on the mean heat transfer. We employ the integral method to obtain the thermal 
solution and utilize the previously derived velocity profiles. 


Integral Relations 


The boundary -layer heat transfer relies mainly on the steady temperature distribu- 
tion due to the basic mean flow convection at constant plate -to -stream temperature dif- 
ference. Temperature oscillations are of the first order (in e) and are caused by the 
convective effect of the first-order velocity fluctuations. The steady second-order tem- 
perature is a result of the interaction between velocity and temperature oscillations, 
plus the convective effect of the steady mean flow correction u 2 - With these consider- 
ations, the dimensionless temperature is expanded as follows: 


6 = 0 Q (x,y) + e^ 1 (x,y)e iajt + 


e OJx, y) + 0«(x,y,t) 


+ ^(e 3 ) 


(75) 


We substitute equation (75) into the thermal integral relation (14), nondimension- 
alize, and use the velocity profile expansion (16) to obtain the following thermal integral 
relations: 

Order : 
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Order 
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( 78 ) 



where Pr = v/oi^ is the Prandtl number. 

We again use limiting -case solutions as guidelines to assume approximate 
temperature -profile forms. The detailed temperature distributions are then deter- 
mined to satisfy the thermal integral equations (76) to (78) . 


Basic Mean Temperature 


The basic mean temperature distribution is identified with the flat-plate boundary 
layer in a uniform free stream with zero plate temperature and unit stream tempera- 
ture. This temperature distribution is the well-known similarity solution of the follow- 
ing boundary -value problem (cf. ref. 12): 


G" + iPrF^G’ - 0 
2 B 

G(0) = 0 


G(tj) -1 as ))-“ J 


(79) 


The resulting solution is 


l [ F B< z >] Pr,iZ 

G(v) = (80) 

/> WF dz 

J o 


We obtain an approximation for the profile (80) by assuming the temperature dis- 
tributions of the following form: 


6 q = 1 - e~ Ar] (l + Bt? + Crj 2 + Dr? 3 ) 


(81) 


The constants A, B, C, and D are computed to satisfy the thermal integral rela- 
tions (76) -with k = 0, 1, 2, and 3. Upon changing variables ( 77 = and introducing 

the function notation, there results 


[F(3, 2) -F(3,l)] + A (a - b)= 0 
Pr ' 


( 82 ) 
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(83) 


[F(3, 3) - F(3, 3, 2) - 2F(3, 1) + 2F(3, 2)] + -± [F(4, 4)] = 0 

Pr 

[F(3, 3, 3, 2) - F(3, 3, 3) - 3F(3, 3, 2) + F(3, 3) + F(3, 2) - 3F(3,1)] 

+ — [F(4, 4) - F(4, 4, 3)] = 0 (84) 

Pr 

[F(3, 3, 3, 3) - F(3, 3 , 3, 3, 2) + 4F(3, 3, 3, 2) - 4F(3, 3, 3) - 6F(3, 3, 2) + 6F(3, 3) + 4F(3, 2) 

- 4F(3, 1)] + — [F(4, 4) - 2F(4,4,3) + F(4, 4, 3, 3)] = 0 (85) 

Pr 

where the following definitions are utilized: 
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~0~ 
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0 ~ 

1 - U A 

[ s i] = 
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[Vj - 
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d0 o /d?j 
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A - B 

AB - 2C 

AC - 3D 

AD 


The solution to the set of algebraic equations (82) to (85) is, for Pr = 0.72 


A 


1.665 

B 


1.367 





C 


0.6699 





D 


0.7161 


(87) 


The corresponding profile results are compared with the exact similarity solution in 
table II and also shown graphically in figure 18. In addition to the obvious agreement, 
we note the following comparison of slopes at the wall: 




= 0.2976 1 


> 


= 0.2980 


(88) 
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TABLE II. - COMPARISON OF ZEROTH-ORDER 


TEMPERATURE PROFILE APPROXIMATION 
WITH EXACT QUASI -STEADY SOLUTION 


Similarity 

variable, 

i> = y/S 

Exact 

solution 

Approximate 

profile 

Percent 

error 

0.20 

0.05913 

0.06395 

8.15 

.40 

.11822 

.12670 

7.18 

.60 

.17720 

.18395 

3.81 

.80 

.23590 

.23746 

.66 

1.00 

.29416 

.28997 

-1.43 

1.20 

0.35171 

0.34333 

-2.38 

1.40 

.40828 

.39817 

-2.48 

1.60 

.46355 

.45412 

-2.03 

1.80 

.51716 

.51025 

-1.34 

2.00 

. 56876 

.56539 

-.59 

2.20 

0.61800 

0.61838 

0.06 

2.40 

.66455 

.66828 

.56 

2.60 

.70812 

.71439 

.89 

2.80 

.74845 

.75626 

1.04 

3.00 

.78536 

.79370 

1.06 

3.20 

0.81875 

0.82672 

0.97 

3.40 

.84858 

.85549 

.81 

3.60 

. 87490 

.88027 

.61 

3.80 

.89781 

.90140 

.40 

4.00 

.91750 

.91927 

.19 

4.20 

0.93418 

0.93424 

0.01 

4.40 

.94812 

.94670 

-.15 

4.60 

.95962 

.95699 

-.27 

4.80 

.96896 

.96544 

-.36 

5.00 

.07645 

.97235 

-.42 

5. 20 

0.98236 

0.97795 

-0.45 

5.40 

.98697 

.98248 

-.45 

5.60 

.99049 

.98613 

-.44 

5.80 

.99315 

.98905 

-.41 

6.00 

.99515 

.99138 

-.38 

6.20 

0.99662 

0.99323 

-0.34 

6.40 

.99769 

.99470 

-.30 

6.60 

.99845 

.99586 

-.26 

6.80 

.99899 

.99678 

-.22 

7.00 

.99936 

.99749 

-.19 





Figure 18. - Four-parameter profile approximation of the 
basic mean temperature distribution. 


The wall slope is, of course, related to the basic mean heat transfer to the wall. 


Temperature Oscillations 

Oscillations in the temperature distribution are generated by the oscillating compo- 
nent of velocity. We compute these fluctuations by applying the integral method of 
analysis. 

Quasi- steady temperature . - We first consider the quasi- steady effect of low- 
frequency velocity fluctuations, that is, the steady temperature distribution with in- 
stantaneous uniform stream velocity Uq( 1 + ee iwt ), namely, 



Now, the temperature distribution (89) is expanded in a Taylor series for small free- 
stream perturbations and the first-order quasi-steady temperature is readily identified 
as the following : 

®1Q = — (90) 

2 

This profile (for Pr = 0.72) is indicated by the dashed line in figure 19. 



Figure 19. - First-order quasi-steady temperature distribution. 


An approximation for is obtained with the integral method by assuming the pro- 
file form 


= e -aT? (b ?7 + cr^ + d^) 


(91) 


39 


We introduce the following definitions: 
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— _ 











Then, following the usual procedure, the thermal integral relations (77) for k = 0, 1, 
2, and 3 are reduced, for the quasi-steady case, to the following : 


(92) 


[F(6, 1) - F(6, 2)] + [F(5, 3) - F(3, D] + = 0 (93) 

Pr 

2[F(6, 1) - F(6, 2) - F(6, 3) + F(6, 3, 2)] + [F(3, 3) - F(5, 3, 3) - 2F(3, 1) 

+ 2F(5,3)] + A[F(7,4)] = 0 (94) 

Pr 

3[F(6, 1) - 2F(6, 3) - F(6, 2) + F(6, 3, 3) + 2F(6, 3, 2) - F(6, 3, 3, 2)] + H [F(6, 4, 4) 

Pr 

+ 2F(7, 4) - 2F(7, 4, 3)] - [F(3, 3, 3) - 3F(3, 3) + 3F(3, 1) 
- F(5, 3, 3, 3) + 3F(5, 3, 3) - 3F(5, 3)] - 0 (95) 

4[F(6, 1) - 3F(6, 3) + 3F(6, 3, 3) - F(6, 3, 3, 3) - F(6, 2) + 3F(6, 3, 2) - 3F(6, 3, 3, 2) 

+ F(6, 3, 3, 3, 2)] + [F(3, 3, 3, 3) - 4F(3, 3, 3) + 6F(3, 3) - 4F(3, 1) - F(5, 3, 3, 3, 3) 

+ 4F(5, 3, 3, 3) - 6F(5, 3, 3) + 4F(5, 3)] + — [F(7, 4) - 2F(7, 4, 3) + F(7, 4, 3, 3) 

Pr 

+ F(6,4,4) -F(6, 4,4,3)] = 0 (96) 
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The solution of the profile constants from these equations gives 


a 


1. 31812 

"b 


0.14310 

"c 


-0.11436 

"d 


0. 50646_ 


( 97 ) 


The resulting quasi-steady profile approximation is shown in figure 19. We also note 
here the corresponding wall slopes for the exact limiting -case solution and approximate 
result: 





0.1488 


> 



0.1431 


(98) 


Profile parameters . - The general fluctuating temperature profile is assumed to be 
of the following form: 

flj. = e~ ay (& + yy 2 + 6y 3 ) (99) 

The coefficients (a, /3, y, and 5) are complex functions of the variable C = ywX/Ug. 
Thus, the temperature oscillations have components in phase and out of phase with the 
free-stream velocity fluctuations. 

The thermal integral relations (77) reduce to a set of ordinary differential equations 
when the profile (99) is introduced. The details of the differential equations are given in 
appendix E and the equations are of the form 


P il- + P i2^ +P i3^ + P i4- =R i 
11 d? ^d? u d? 14 d£ 1 


i = 1, 2, 3, 4 


( 100 ) 


Of course, the real initial values for equation (100) are the quasi-steady limits, namely, 
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as C — 0 


( 101 ) 
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_ r_ 


dA 3 


We obtain the imaginary initial values by introducing an expansion of the identical form 
as equations (59). Substituting the expansions into equation (100) and solving the result- 
ing algebraic equations by the method of appendix D then yields 


“i 


0.06196 £ 

A 


-0.02901 C 

A 


-0.08295 

6. 

lj 


0.0688 5 A 


(102) 



Figure 20. - Real parts of first-order temperature profile form parameters. 
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Figure 21. - Imaginary parts of first-order temperature profile form parameters. 


Numerical solution of equations (100) with initial values (101) and (102) gives the 
profile form parameters shown graphically in figures 20 and 21 . 

Temperature profile results . - The in-phase and out -of -phase temperature profiles 

are plotted against j] in figures 22 and 23 for small values of £ = -^/wX/Up. We see 
that as the profiles develop, the peak value increases but remains at approximately the 
same value of tj. 

The profile development for large values of £ is shown in figures 24 and 25 with the 
independent variable y. We note the spreading of the effected layer into the free stream 
and ultimately a reduction of both components for very large £ . 

We define the engineering heat -transfer quantities in the usual manner. The heat 
flux at the wall is 




• k (-L 


The flux is related to the Nusselt number Nu^ as follows: 


Qw = | Nu X< T w- T ^ 


(103) 


(104) 
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We consider here the ratio of the first-order to the mean-flow Nusselt number and de- 
fine this heat -flux ratio as follows: 



(105) 


The real and imaginary parts of the heat-flux ratio are shown in figure 26. 

The first-order heat transfer is oscillatory in nature and thus gives no net contri- 
bution to the overall heat transfer in the time -average sense. 




(b) Imaginary part. 

Figure 26. - First-order wall heat-flux ratio. 
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Effect of Oscillations on the Time-Average Temperature 

We now compute the effect of velocity and temperature oscillations on the time- 
average temperature distribution. The steady mean temperature correction is a result 
of the interaction between the first-order temperature and velocity oscillations plus the 
thermal convection due to the second -order mean flow. 

Quasi-steady temperature . - The second-order quasi-steady temperature distribu- 
tion (limiting case for low-frequency fluctuations) is computed directly from the small 
perturbation expansion of equation (89). The result is 


«2Q = - O’) 


(106) 


This profile is given in figure 27. 



Figure 27. - Second-order quasi-steady temperature profile. 
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The quasi-steady approximation profile is assumed to have the following form: 


0 2 q = e" pT? (qrj + rr/ 2 + st] 2 ) 


(107) 


We now introduce the arrays 
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Then, following the usual procedure for similarity profiles, we reduce the thermal inte- 
gral relations (78) for k = 0, 1, 2, and 3 to the following: 


[F(6, 1) - F(6, 3)] + 2[F(7, 1) - F(7, 2)] - 2[F(5, 4)] + -i- -q = 0 (109) 

Pr 

[F(6, 6) - F(6, 6,2)] +4[F(7, 1) - F(7, 2) - F(7, 5) + F(7, 5, 2)] + 2[F(6, 1) - F(6, 5) - F(6, 3) 

+ F(6, 5, 3)] + 2[F(5, 5, 4) - 2F(5, 4)] + A [ 4 F(10, 8) + F(9, 9)] - 0 (110) 

Pr 
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3[F(6,6) - F(6, 6, 2) - F(6, 6, 5) + F(6, 6, 5, 2)] + 6[F(7, 1) - F(7, 2) - 2F(7, 5) + 2F(7, 5,2) 

+ F(7, 5, 5) - F(7, 5, 5, 2)] + 3[F(6, 1) - 2F(6, 5) - F(6, 3) + 2F(6, 5, 3) + F(6, 5, 5) 

- F(6, 5, 5,3)] - 2[F(5, 5, 5,4) - 3F(5, 5,4) + 3F(5,4)] + 11 [4F(10, 8) 

Pr 

- 4F(10, 8, 5) + 2F(9, 8, 6) + 2F(8,8,7) + F(9, 9) - F(9,9,5)] = 0 (111) 

3[F(6, 6) - F(6, 6, 2) - 2F(6, 6, 5) + 2F(6, 6, 5, 2) + F(6, 6, 5, 5) - F(6, 6, 5, 5, 2)] + 4[F(7, 1) 

-F(7, 2) - 3F(7, 5) + 3F(7, 5, 2) + 3F(7, 5, 5) - 3F(7, 5, 5, 2) - F(7, 5, 5, 5) 

+ F(7, 5, 5, 5, 2)] + 2[F(6, 1) - 3F(6, 5) - F(6, 3) + 3F(6, 5, 3) + 3F(6, 5, 5) 

- 3F(6, 5, 5, 3) - F(6, 5, 5, 5) + F(6, 5, 5, 5, 3)] - [4F(5, 4) - 6F(5, 5, 4) 

+ 4F(5, 5, 5, 4) - F(5, 5, 5, 5, 4)] + H [4F(10, 8) - 8F(10, 8, 5) + 4F(10, 8, 5, 5) 

Pr 

+ 4F(9, 8, 6) - 4F(9, 8, 6, 5) + 4F(8, 8, 7) - 4F(8, 8, 7, 5) + F(9, 9) - 2F(9, 9, 5) 

+ F(9,9, 5,5) + F(8, 8,6,6)] = 0 (112) 


The preceding equations admit the following solution: 


p 


" 1.08776 " 

q 


-0.018378 

r 


0.091511 

s 


_-0. 12387 


(113) 


The profile corresponding to the result (113) is shown in figure 27. The approximation 
gives the general shape desired, and the peak values agree very well in magnitude and 
reasonably well in location. 

It is of interest that the slope of the approximate profile at the wall compares with 
the quasi-steady limit value as follows: 
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= -0.01860 






> 


-0.01838 


j 


(114) 


These slopes agree even though the curves have different shapes near the wall. We in- 
vestigate this effect by assuming profiles with other slopes at the wall (i.e. , pick q in 
eq. (107)). The remaining three shape factors in equation (107) are determined to 
satisfy equations (110) to (112). The profiles considered are 


0 2 ^ = e -1 ' 076lT7 (-0.10 rj + 0.07370 tj 2 
0^ = e - 1,0639t? (-0.30 n + 0.07024 tj 2 


A 


0.1154 v ) 


0.1084 ^ 


(115) 


The behavior of the various profiles is compared in figure 28. The profile value for 
rj > 1 is completely determined by the coefficients p, r, and s. Thus, the value of q 
which satisfies equation (109) is that which causes 6 ^ to satisfy the heat -flux integral 
near the wall. 

Profile factors . - The steady second-order temperature profile, which applies for 
all frequency values, has the form 


6 2 = e -<iPy (^y + jzy 2 + Xy 3 ) 


(116) 


We compute the profile form factors in the usual manner by the integral method. The 
thermal integral relations (78) reduce to differential equations of the following form: 



~ +P i2 — + P i3^ L+ P i4~ = R i 
d? * d? u d? 14 d? 1 


i= 1, 2, 3, 4 


(117) 


The coefficients P^ and 


are listed in appendix E . 
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The initial values of the profile factors are obtained from the quasi-steady distribu- 
tion, that is, 



The numerical solution of this initial -value problem, equations (117) and (118), is pre- 
sented graphically in figure 29 . 

Temperature distribution . - The development of the steady second-order tempera- 
ture profiles is shown in figures 30 and 31. We see that the initial departure from 



Figure 30. - Steady second-order temperature profiles for F '9 ure 31 - ' stead y second-order temperature profiles for large values 

small values of C of ?• 
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Result of - 


Prandtl 

number 


Present study 0.72 

Reference 8 . 70 

Reference 7 (+Cq 2 .70 

shown including 
effects of U2 



Figure 32. - Steady second-order heat-flux ratio. 


quasi-steady behavior is to increase the maximum value, with that maximum occurring 
at about the same value of rj (fig. 30). The layer then spreads into the stream (fig. 31) 
and ultimately diminishes in value. 

A second-order heat -flux ratio is defined as follows: 



(119) 


Values of this ratio are compared with results of other analytical investigations in fig- 
ure 32. All previous work is asymptotically valid for small o>X/Uq . Thus, the trend 
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of the Gersten result is confirmed over the full range of frequency. The Nickerson data 
are for the heat -transfer contribution of the second-order steady velocity and neglect 
the oscillating velocity -temperature interaction. 

We use the preceding heat -transfer results, together with equation (103), to express 
the mean rate of heat transfer to the plate. The result for Pr = 0. 72 is 



( 120 ) 


The corresponding values of Cqg are shown in figure 32. Since Cqg is negative, the 
effect of flow oscillations is a small reduction in the mean heat transfer to the plate. 


SUMMARY AND CONCLUSIONS 

The general integral method of analysis was developed for the unsteady thermal 
boundary layer. The development was based on the flat -plate geometry with time- 
dependent free-stream velocity. The stream and plate were taken to be at constant but 
different temperature levels. 

A system of integral relations was derived based on the momentum and thermal 
energy equations. The flow integral relations are velocity -weighted averages of the 
momentum equation. The thermal integral relations are temperature -weighted averages 
of the thermal energy equation. 

The integral method then consists of the following: Obtain asymptotic solutions for 
limiting cases of the flow and heat transfer under consideration. Assume velocity and 
temperature profile forms with sufficient generality to include the asymptotic solutions. 
Determine profile form factors to satisfy the appropriate integral relations. 

A detailed analysis was performed for the special case of a free stream with small 
harmonic velocity oscillations about a steady mean value. The velocity was found to 
consist of a steady mean flow corresponding to the Blasius flow, first-order velocity os- 
cillations, and a second -order correction to the mean flow due to feedback of the oscil- 
lations. The basic mean temperature distribution was found to be the well-known flat- 
plate steady -boundary -layer result. First-order temperature oscillations were gener- 
ated by the velocity oscillations. 

The second-order steady temperature correction was caused by the velocity- 
temperature oscillation interaction plus the steady second-order flow convection effect. 
Oscillations were found to reduce the mean heat transfer to the plate. The amount of 
decrease is represented by a single coefficient, which is presented graphically. Flow 


54 



and heat-transfer results were found to adequately span the asymptotic limiting -case 
solutions. 

Although the analysis herein was specialized for certain geometry, free stream, 
and boundary conditions, the basic method could be applied to many other configurations 
of interest. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, August 4, 1972, 

503-10. 
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SYMBOLS 

Blasius profile approximation form factor 
zeroth-order temperature profile form factor 
first-order quasi-steady velocity profile factor 
first-order quasi-steady temperature profile factor 
imaginary coefficient in a expansion, eq. (59) 
Blasius profile approximation form factor 
zeroth-order temperature profile form factor 
first-order quasi -steady velocity profile factor 
first-order quasi-steady temperature profile factor 
imaginary coefficient in (3 expansion, eq. (59) 
Blasius profile approximation form factor 
zeroth-order temperature profile form factor 
constant, eq. (36) 
heat-flux ratio 

first-order wall shear stress coefficient, eq. (62) 
second -order wall shear stress coefficient, eq. (74) 
first-order quasi-steady velocity profile factor 
first-order quasi-steady temperature profile factor 
imaginary coefficient in y expansion, eq. (59) 
Blasius profile approximation form factor 
zeroth-order temperature profile form factor 
first-order quasi-steady velocity profile factor 
first-order quasi-steady temperature profile factor 
imaginary coefficient in 5 expansion, eq. (59) 
Blasius flow dimensionless stream function 
profile summation function 

zeroth-order flat -plate temperature profile, eq. (80) 
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lm 
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1m 
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Nu x 


n i 



P 

P 

Q 

q 

q 

R. 

1 

r 

r 


s 

s 


T 

t 

U 

u o 

u 


functions, appendix D 
profile approximation matrix 
profile approximation matrix 
integrals, appendixes B and C 
nonnegative integer 
Runge-Kutta functions, appendix F 
Runge-Kutta functions, appendix F 
Runge-Kutta functions, appendix F 
Nusselt number 

Runge-Kutta functions, appendix F 
pressure 

coefficients in differential equations 
Prandtl number , 

second -order quasi -steady velocity profile factor 
second -order quasi-steady temperature profile factor 
heat flux at wall 

second-order quasi-steady velocity profile factor 
second-order quasi-steady temperature profile factor 
terms in differential equations 
second-order quasi-steady velocity profile factor 
second -order quasi-steady temperature profile factor 
attenuation matrix 

second-order quasi-steady velocity profile factor 

second-order quasi-steady temperature profile factor 

attenuation matrix 

temperature 

time 

velocity component parallel to the body 
mean value of free -stream velocity 

dimensionless velocity component parallel to the body, U/Uq 
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u 2 

(N» 

u 2 

V 

X 

x 

Y 

y 

z 

a 

ct 

“th 

/3 

£ 

y 

y 

A 

6 

6 

€ 

s 

77 

e 

K 

X 

x 

fj- 

M 

y 

7T 


time -independent part of second -order velocity 
time -dependent part of second-order velocity 
velocity component normal to the body 
coordinate along the body 

dimensionless coordinate along the body, wX/Uq 
coordinate normal to the body 


dimensionless coordinate normal to the body, 



Y 


dummy variable, eq. (36) 
first-order velocity profile form function 
first-order temperature profile form function 
thermal diffusivity 

first-order velocity profile form function 
first-order temperature profile form function 
first-order velocity profile form function 
first-order temperature profile form function 
increment in Runge-Kutta integration formulas, appendix F 
first-order velocity profile form function 
first-order temperature profile form function 
amplitude of harmonic free -stream velocity fluctuation 
stretched coordinate, V* 
similarity variable, y/£ 

dimensionless temperature, (T - TJ/(T„ - T ) 

thermal conductivity 

second -order velocity profile factor 

second-order temperature profile factor 

second -order velocity profile factor 

second -order temperature profile factor 

kinematic viscosity 

differential function, eq. (57) 
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p 

density 

*r 

phase angle of first-order wall shear stress coefficient 

9 

second-order velocity profile form function 

<P 

second -order temperature profile form function 


second -order velocity profile form functions 

V 

second -order temperature profile form functions 

S2 

function, eq. (25) 

u > 

frequency of harmonic free -stream fluctuations 

Subscripts: 


B 

Blasius flow value 

i 

imaginary part 

Q 

quasi -steady expansion value 

q 

quasi-steady approximation value 

r 

real part 

S 

Stokes flow value 

w 

value at plate surface 

0 

zeroth-order (e®) value 

1 

first-order (e 1 ) value 

2 

second -order (e 2 ) value 

OO 

free -stream value 

Special symbols: 

f(i,j,k,l) 

integration summation function, eq. (35) 

G(I,J,K,L) 

integration summation function, eq. (E2) 

C 

order symbol 

Re{ } 

real part 

Im{ } 

imaginary part 

Z* 

denotes complex conjugate of Z, or prescribed Z in appendix D 

|z| 

denotes magnitude of Z 

g<+) 7(") 

temperature profiles, eq. (115) 
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APPENDIX B 


DERIVATION OF THE FLOW INTEGRAL RELATIONS 

The X-momentum equation (2) is multiplied by a nonnegative power of the local axial 
velocity U k and integrated in the direction normal to the plate from the surface to the 
free stream. The result is 



P 00 

/ " 

Jo 


k 

at 


dY + v 


I, 


U k ^dY 


(Bl) 


0 ZY‘ 


First, consider the following integral: 


I 


1 , = I VU k — dY 

1 1 dY 


(B2) 


We note that the V-velocity is related to the U-velocity by the continuity equation (1). 
Integrating the continuity relation in the normal direction gives 



Then equation (B2) can be written as follows: 
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Integrating equation (B4) by parts we obtain 



Next consider the integral 




dY 


For k = 0, the integral is evaluated directly and the result is 


(B5) 


(B6) 


M° = .(3U\ k = 0 

\3Y/ y= 


For k ^ 0, equation (B6) is integrated by parts and we obtain 


(Bl) 



Finally, upon substituting the results (B5), (B7), and (B8) into equation (Bl) , we ob- 
tain the following general integral relation: 
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APPENDIX C 


DERIVATION OF THE THERMAL INTEGRAL RELATIONS 


We multiply the thermal energy equation (4) by a nonnegative power of the dimen- 
sionless temperature and integrate the result in the Y-direction from the plate surface 
to the free stream. The result is 



The transverse velocity V is eliminated by applying the continuity requirement, namely, 



Then equation (C2) becomes 



Integration of equation (C4) by parts yields 
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Next consider the following: 



For k = 0, direct evaluation of the integral gives the result 


I 


4 


( d9 \ k 

W Y= 


(C6) 


(C7) 


For k / 0, we integrate equation (C6) by parts and the result is as follows: 



Finally, upon substituting equations (C5), (Cl), and (C8) into the integrated energy 
equation (Cl), we obtain the following thermal integral relations: 
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APPENDIX D 


METHOD OF SOLUTION FOR SYSTEMS OF NONLINEAR ALGEBRAIC EQUATIONS 

The determination of the velocity profile form parameters often requires solving a 
system of nonlinear algebraic equations of the following form: 

gj(A, B,C, . . .) = 0 (Dl) 

In this appendix, the method employed in this work for all such calculations is presented. 
First, consider the simple case of one equation of the type (Dl); that is, 


gj(A) = 0 


(D2) 


In this case, a graph of gj against A is constructed and the appropriate root is 
selected . 

Next consider two equations of the following form: 


gj(A , B) 
g 2 (A,B) 


= 0 

(D3) 

= 0 

(D4) 


A value of A (A = A*) is picked, and the first equation is evaluated for that A* for two 
selected values of B; that is, 




(D5) 

(D6) 


Referring to figure 33, the value of B for which gj is zero (with A = A*) is computed 
by using the following approximation: 



(D7) 
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Figure 33. - Linear approximation for a nonlinear function of one variable. 


The value of gt,(A* ,B*) is computed and graphed against A* to determine the solution 
set (A,B). The process is repeated, with selected values of B chosen closer to the 
calculated root until the desired accuracy is achieved. 

A system of three equations of the type (Dl) has the form 


gl (A,B,C) = 0 

(D8) 

g 2 (A,B,C) = 0 

(D9) 

g 3 (A,B,C) =0 

(DIO) 


A procedure is employed which is an extension of the above procedure for two equations. 
Assumed values of C are also required, and the linear representations shown in fig- 
ure 34 are utilized. The resulting equations for B* and C* are the following: 


C* 



B* 




(Dll) 


(D12) 


65 






Now, a graph of gg against A* is constructed to obtain the roots, and the procedure 
is continued to achieve the desired accuracy. 

The procedure is directly extended for larger values of i in equation (Dl). 
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APPENDIX E 


DIFFERENTIAL EQUATIONS FOR PROFILE PARAMETERS 

In this appendix, we list coefficients of the differential equations of the following 
form: 



(El) 


Values corresponding to the first- and second-order flow and temperature approxima- 
tions are given. All cases utilize the following definitions: 




H t (l + m - 2)'. 

G(L,M 


l sum on l , m = 1 , 2, . . . , 5 


(E2) 


( S L + S n) 


Z+m-1 


First-Order Velocity 

For the first-order oscillating velocity profiles, the following definitions apply: 
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1 


" 0 


1 

0 

0 

0 

0 

1 - u 0 


A/C 


1 

B/5 

c/5 2 

D/5 3 

0 

3u Q /3y 


A/C 


(A - B)/r 

(AB - 2C)/f 2 

(AC - 3D)/? 3 

AD/5 4 

0 

1 - u i 


a 


1 

0 

y 

6 

0 

3 Uj/3 y 


a 


a - 0 

a/3 - 2 y 

ay - 36 

aid 

0 

A J 

[S t ] = 




-1 

-(3 


-6 

-3Uj/3a 

Ct 

0 

-y 

-BUj/3/3 


a 


0 

1 

0 

0 

0 

-Suj/ay 


a 


0 

0 

1 

0 

0 

-3Uj/36 


a 


0 

0 

0 

1 

0 

-3u 0 M 


A/C 


0 

(A - B)/5 2 

(AB - 2C)/5 3 

(AC - 3D)/5 4 

AD/5 5 


(E3) 


Then, the coefficients of equation (El) which are utilized in equation (55) are 

P lj = G(j + 5,l) - 2G(j + 5,2) 

P 2j = 2G(j + 5,1)- 6G(j + 5, 2) + 3G(j + 5, 2, 2) 

P 3j = 3G(j + 5,1) - 12G(j + 5,2) + 12G(j + 5,2,2) - 4G(j + 5, 2, 2, 2) 
P 4j = 4G(j + 5,1)- 20G(j + 5,2) + 30G(j + 5, 2, 2) - 20G(j + 5, 2, 2, 2) 

+ 5G(j + 5, 2, 2, 2, 2) 




> 




(E4) 


R x = 2G(10, 4) - 2?iG(4,l) + 2?(a - j3) - F(2,l) (E5) 

R 2 = 6[G(10, 4) - G(10, 4, 2)] - 4fi[G(4, 1) - G(4, 2)] + 8CG(5, 3) - [4F(2, 1) - 3F(2, 2)] (E6) 

R 3 = 12[G(10, 4) - 2G(10, 4, 2) + G(10, 4, 2, 2)] - 6?i[G(4, 1) - 2G(4, 2) + G(4, 2, 2) ] 

+ 12?[2G(5, 3) - 2G(5, 3, 2) - G(4, 3, 3)] - [9F(2, 1) 

- 12F(2, 2) + 4F(2, 2, 2) - 12? 2 F(3, 3)] (E7) 
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R 4 = 20[G(10, 4) - 3G(10, 4, 2) + 3G(10, 4, 2, 2) - G(10, 4, 2, 2, 2)] - 8?i[G(4, 1) j- 3G(4, 2) 

+ 3G(4, 2, 2) - G(4, 2, 2, 2)] + 48?[G(5, 3) - 2G(5, 3, 2) + G(5, 3, 2, 2) 

- G(4, 3, 3) + G(4, 3, 3, 2)] - [16F(2, 1) - 30F(2, 2) + 20F(2, 2, 2) 
-5F(2, 2, 2, 2) - 48£ 2 G(3, 3) + 48? 2 G(3, 3, 2)] (E8) 

The F -function definitions corresponding to equation (48) have been retained. 


Time-Independent Second-Order Velocity 

For the time -independent second -order velocity profile approximations the follow- 
ing arrays are utilized: 
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0 

Q 

0 
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1 -Uq 
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1 
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CA 2 

DA 3 

0 

1-Uj 
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1 

0 

r 

6 

0 
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a * 



1 

0* 

7* 

6* 

0 
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0 

* 

M 

A 

0 

3u 0 ,/3y 


A A 


(A 

- B)A 

(AS - 2C)/< 2 
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ADA 4 

0 

3u 1 /3y 
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ay - 36 

ff6 

0 

3uJ /3y 


a* 
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a* 0* - 2y * 

a*y* - 36* 
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0 
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<P 
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0 
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0 
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Then, the coefficients of equation (73) have the values 




V 1 . =[G(j + 12, 1) - 2G(j + 12, 2)] 

P 9 . = I [2G(j + 12, 1) - 6G(j + 12, 2) + 3G(j + 12, 2, 2)] 

^ 2 

P 3 . = I[3G(j + 12, 1) - 1 2G(j + 1 2, 2) + 1 2G(j + 1 2, 2, 2) - 4G(j + 1 2, 2, 2, 2) ] 
3 

P 4 . = I [4G(j + 12, 1) - 20G(j + 12, 2) + 30G(j + 12, 2, 2) 

4 

- 20G(j + 12, 2, 2, 2) + 5G(j + 12, 2, 2, 2, 2)] 


> (E10) 


R 1 = - 2G(10, 5) + - Re {G(ll, 4) - G(12, 1) + G(12, 3)} 

1 2 


(Ell) 


R 2 = -3[G(10, 5) - G(10, 5, 2)] - SJ?{G(3, 1)} - 4£G(9, 6) - $ Re {G(8, 7)} - - Re {2G(11, 1) 

4 

+ 2G(10, 1) + 3[G(12, 1) - G(12, 3) - G(12, 2) + G(12, 3, 2) - G(ll, 4) 
- G(ll, 2) + G(ll, 4, 2) - G(10, 4) - G(10, 3) + G(10, 4, 3)]} (El 2) 
R 3 = -4[G(10, 5) - 2G(10, 5, 2) + G(10, 5, 2, 2)] - 2^{G(3, 1) - G(3, 2)} - 8?[G(9, 6) 

- G(9, 6, 2)] - 2C &*{G(8,7) - G(8, 7, 2) + G(8, 6) - G(8, 6, 3) + G(7, 6) - G(7,6,4) 

+ 2G(6, 6, 5)} - Re G(10, 1) - 2G(10, 2) - 2G(10, 3) - 2G(10, 4) + 2G(10, 4, 3) 

+ 2G(10, 4, 2) + 2G(10, 3, 2) - 2G(10, 4, 3, 2) + G(12, 1) - G(12, 3) - 2G(12, 2) 

+ 2G(12, 3, 2) + G(12, 2, 2) - G(12, 3, 2, 2) + ~ G(ll, 1) - G(ll, 4) - 2G(11, 2) 

2 

+ 2G(11, 4, 2) + G(ll, 2, 2) - G(ll, 4, 2, 2)J (El 3) 
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R 4 = -5[G(10, 5) - 3G(10, 5, 2) + 3G(10, 5, 2, 2) - G(10, 5, 2, 2, 2)] - 3^{G(3, 1) - 2G(3, 2) 

+ G{3, 2, 2)} - 12?[G(9, 6) - 2G(9, 6, 2) + G(9, 6, 2, 2) + G(6, 6, 5) - G(6, 6, 5, 2)] 

- 3? Re { [G(8 , 7) - 2G(8, 7, 2) + G(8, 7, 2, 2)] + 2[G(8, 6) - G(8, 6, 2) 

-G(8,6,3) + G(8,6,3,2)]+ 2[G(7, 6) - G(7, 6, 2) - G(7, 6, 4) + G(7, 6, 4, 2)] 

+ [G(6,6) - G(6, 6, 4) -G(6,6,3) + G(6, 6, 4, 3)]}- - Re j3[G(10, 1) 

- 2G(10, 2) + G(10, 2, 2) - G(10, 3) + 2G(10, 3, 2) - G(10, 3, 2, 2) - G(10, 4) 

+ 2G(10, 4, 2) - G(10, 4, 2, 2) + G(10, 4, 3) - G(10, 4, 3, 2) + G(10, 4, 3, 2, 2)] 

+ G(12, 1) - 3G(12, 2) + 3G(12, 2, 2) - G(12, 2, 2, 2) - G(12, 3) + 3G(12, 3, 2) 

- 3G(12, 3, 2 , 2) + G(12, 3, 2, 2, 2) - ^ G(ll, 1) - 3G(11, 2) + 3G(11, 2, 2) 

5 

- G(U, 2, 2, 2) - G(ll, 3) + 3G(11, 3, 2) - 3G(11, 3, 2, 2) + G(ll, 3, 2, 2, 2) J 

(El 4) 


First-Order Temperature 

The first-order temperature profile form factor calculation is based on the follow- 
ing definitions: 
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With these definitions, the terms in equation (100) are 


P Xj =[G(j + 12,1) - G(j + 12,2)] 

P 2j = 2[G(j + 12,1) - G(j + 12, 2) - G(j + 12,8) +G(j + 12,8,2)] 
P 3j = 3[G(j + 12, 1) - G(j + 12, 2) - 2G(j + 12, 8) + 2G(j + 12, 8, 2) 

+ G(j + 12,8,8) - G(j + 12,8,8,2)] 
P 4j = 4[G(j + 12, 1) - G(j + 12, 2) - 3G(j + 12, 8) + 3G(j + 12, 8, 2) 
+ 3G(j + 12, 8, 8) - 3G(j + 12, 8, 8, 2) - G(j + 12, 8, 8, 8) 

+ G(j + 12,8,8,8,2)] 




J 


(E15) 


(El 6) 
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4 = - - 2?i[G(ll, 1)] - [G(ll, 4) + G(10, 1) - G(10, 5) - G(8, 7)] (El 7) 


R 2 = ~(^") [G(12 ’ 9)] " 4 ^ G(11 > 1} ' G(11 ’ 8) J ‘ 2 f G(11 ’ 10) " G(11 ' 10 ’ 2) + G(U > 4 > 


- G(ll, 8, 4) + G(10, 1) - G(10, 8) - G(10, 5) 


R 3 = [2G(12, 9) - 2G(12, 9, 8) + G(ll,9,9)] - 6?i[G(ll, 1) - 2G(11,8) + G(ll,8,8)] 

- 6[G(11, 10) - G(ll, 10, 2) - G(ll, 10, 8) + G(ll, 10, 8, 2)] - 3[G(11, 4) + G(10, 1) 

- G(10, 5) - 2G(11, 8, 4) - 2G(10, 8) + 2G(10, 8, 5) + G(ll, 8, 8, 4) + G(10, 8, 8) 


+ G(10, 8, 5)] + [G(8, 8, 7) - 2G(8, 7)] (E18) 


(E19) 


- G(10, 8, 8, 5)] - [G(8, 8, 8, 7) - 3G(8, 8, 7) + 3G(8, 7)] 

R,= _[£§£] [ G ( 12>9) - 2G(12, 9, 8) + G(12, 9, 8, 8) +G(11,9,9) - G(ll, 9, 9, 8)] 
4 \Pr/ 


- 8?i[G(ll,l) - 3G(11, 8) + 3G(11, 8, 8) - G(ll,8,8,8)] - 12[G(11,10) - 2G(11,10,8) 
+ G(ll, 10, 8, 8) - G(ll, 10, 2) + 2G(11, 10, 8, 2) - G(ll, 10, 8, 8, 2)] - 4[G(10, 1) 

- 3G(10, 8) + 3G(10, 8, 8) - G(10, 8, 8, 8) - G(10, 5) + 3G(10, 8, 5) - 3G(10, 8, 8, 5) 

+ G(10, 8, 8, 8, 5) - G(ll, 4) - 3G(11, 8, 4) + 3G(11, 8, 8, 4) - G(ll, 8, 8, 8, 4)] 


- [G(8, 8, 8, 8, 7) - 4G(8, 8, 8, 7) + 6G(8, 8, 7) - 4G(8, 7)] 


(E 20) 


Time-Independent Second-Order Temperature 


The required notation for the second -order temperature equation is 
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(E21) 
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Then, the coefficients required for equation (117) are 

P Xj = G(j + 18, 1) - G(j + 18, 2) 

P 2 . = Pj. - G(j + 18, 5) + G(j + 18, 5, 2) 

P 3 . = P^ - 2G(j + 18, 5) + 2G(j + 18, 5, 2) - G(j + 18, 5, 5) 

- G(j + 18, 5, 5, 2) 

P 4j = P Xj - 3G(j + 18, 5) + 3G(j + 18, 5, 2) + 3G(j + 18, 5, 5) 

- 3G(j + 18, 5, 5, 2) -G(j + 18,5,5, 5) 
+ G(j + 18,5,5,5, 2) 

Rj = _ [G( 13 , 8)] - [G(16, 4) - G(15, 5) + 0. 5 G(14, 7) 




>■ 


J 


(E22) 


+ 0.5 G(18, 1) - 0. 5 G(18,3)] (E23) 

R 2 - [G(12, 9) + 0. 25 G(ll, 10)] - [G(13, 8) - G(13, 8, 5) + G(16, 8) - G(16, 8, 2)] 

- [0. 5 G(15, 5, 5) - G(15, 5) + G(16, 4) - G(16, 5, 4)] - 0. 5[G(14, 7) - G(14, 7, 5) 
+ G(18, 1) - G(18, 3) - G(18, 5) + G(18, 5, 3) + G(16, 7) - G(16, 7, 3)] 

- 0. 25[G(13, 7, 6) + G(17, 7) - G(17, 7, 2) + G(18, 6) - G(18, 6, 2)] (E24) 
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R 3 = ^2i^[4G(12,9) - 4G(1 2, 9, 5) + G(11,9,6) + G(10,9,7) + 2G(9, 9, 8) + G(ll,10) 

- G(ll, 10, 5)] - [G(13, 8) - 2G(13, 8, 5) + G(13, 8, 5, 5) + 2G(16, 8) - 2G(16,8, 5) 

- 2G(16, 8, 2) + G(16, 8,5,2)] + — [G(15, 5, 5, 5) - 3G(15, 5, 5) + 3G(15, 5)] 

3 

- [G(16, 4) - 2G(16, 5, 4) + G(16, 5, 5, 4) + G(16, 7) - G(16, 7, 5) - G(16, 7, 3) 

+ G(16, 7, 5, 3) ] - 0. 5[G(14, 7) - 2G(14, 7, 5) + G(14, 7, 5, 5) + G(18, 1) 

- 2G(18, 5) + G(18, 5, 5) - G(18, 3) + 2G(18, 5, 3) - G(18, 5, 5, 3) + G(13, 7, 6) 

- G(13, 7, 6, 5) + G(17, 7) - G(17, 7, 5) - G(17, 7, 2) + G(17, 7, 5, 2) + G(18, 6) 

-G(18,6, 5) - G(18, 6, 2) + G(18, 6, 5, 2) +G(16,7,6) - G(16, 7, 6, 2)] (E25) 

r 4 = ^Lj[4G(12, 9) - 8G(12, 9, 5) + 4G(12, 9, 5,5) + 2G(11,9,6) - 2G(11,9,6, 5) 

+ 2G(10, 9, 7) - 2G(10, 9, 7, 5) + 4G(9, 9, 8) - 4G(9, 9, 8, 5) + G(ll, 10) - 2G(11, 10, 5) 

+ G(ll, 10, 5, 5) + G(9, 9, 7, 6)] - [G(13, 8) - 3G(13, 8, 5) + 3G(13, 8, 5, 5) 

- G(13, 8, 5, 5, 5)] - 3[G(16, 8) - 2G(16, 8, 5) + G(16, 8, 5, 5) - G(16, 8, 2) 

+ 2G(16, 8, 5, 2) - G(10, 8, 5, 5, 2)] - -{[F(15, 5, 5, 5, 5) - 4F(15, 5, 5, 5) + 6F(15, 5, 5) 

4 1 

- 4F(15, 5)] + 4[G(16, 4) - 3G(16, 5,4) + 3G(16, 5,5,4) - G(16, 5, 5, 5, 4)]+ 2[G(14, 7) 

- 3G(14, 7, 5) + 3G(14, 7, 5, 5) - G(14, 7, 5, 5, 5) + G(18, 1) - 3G(18, 5) + 3G(18, 5, 5) 

- G(18, 5, 5, 5) - G(18, 3) + 3G(18, 5, 3) + 3G(18, 5, 5, 3) + G(18, 5, 5, 5, 3)] 

+ 6[G(16, 7) - 2G(16, 7, 5) + G(16, 7, 5, 5) - G(16, 7, 3) + 2G(16, 7, 5, 3) 

- G(16,7, 5,5,3) + G(16, 7, 6) - G(16, 7, 6, 5) - G(16, 7, 6, 2) + G(17, 7, 6, 5, 2)] 

+ 3[G(13, 7, 6) - 2G(13, 7, 6, 5) + G(13, 7, 6, 5, 5) + G(17, 7) - 2G(17, 7, 5) 

+- G(17, 7, 5, 5) - G(17, 7, 2) + 2G(17, 7, 5, 2) - G(17, 7, 5, 5, 2) + G(18, 6) 

- 2G(18, 6, 5) + G(18, 6, 5, 5) - G(18, 6, 2) + G(18,6, 5, 2) - G(18, 6, 5, 5, 2)]| 

(E26) 
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APPENDIX F 


FOURTH-ORDER RUNGE-KUTTA INTEGRATION FORMULAS 

In this appendix we list the Runge-Kutta formulas which were used to solve the sys- 
tems of coupled first-order ordinary differential equations of the form 


— = 0,y, 6, 0 

d? 1 

(FI) 

ir 2 ( ‘ a >P> y > 5 > ® 

d? z 

(F2) 

— = ir 3 (of,^,y,6,0 
d? J 

(F3) 

— = ff 4 (a,0,y, 5, C) 
df[ 4 

(F4) 


The initial conditions are as follows: 
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where 


5 n + l = 6 n + ^ n 0 + 2n l + 2n 2 + n 3 } 


(FI 2) 


k O = Air lK-Mn- 5 n>y 


k l " A n l (“n + 2 k 0’ + 2 Z °’ y ° + 2 m °’ 6n + \ n °’ + \ 

k 2 = A7r l(“n + ^ k l’ Pn + \ l vVn + l m V 6 n + \ n V *n + J A ) 
k 3 = A »lK + k 2’ + Z 2> y n + m 2> 6 n + n 2’ + A) 


(FI 3) 
(FI 4) 

(FI 5) 
(FI 6) 


l 0 ~ A 6 n^n^ 


^1 = A ir 2 rn + 2 k 0’ + g ^0’ y n + g m 0’ 6 n + ^V + \ A ! 


(FI 7) 
(FI 8) 


l 2 A ^2 rn + 2 k l» + ~ *1> r n + 2 m l» 6 n + ~ n l’ + 2 A ) 


(F19) 


Z 3 = A7r 2 (a n + k 2’ ^n W 2’ y n + m 2’ 6 n + n 2’ + A) 


(F20) 


m 0 = Aj 3K>^n> 5 n>y 


m l =A,I 3 a n^V + 2 ^0’ y n + ^ m 0’ 6 n + ^ V + | A 


(F21) 

(F22) 


m 2 = A * 3 (% + ~ k V Pn +1 ~ l V Y n + l m V 6 n + ^ n l> V^ A ) 


(F23) 


m 3 = A7r 3 (a n + k 2’ ^n +i 2’ y n + m 2’ 6 n + n 2’ + A) 


(F24) 
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n 0 A ^^n’^n’^n’ ^n’ V 


<F25) 

"l <*n + j'o- *n + j m 0> 6 n + J n 0' 

K ^) 

(F26) 

n 2 = ’ A, 4(“n t J k r <\i+ J*l> y u + J m l- 5 n + J n l- 


(F27) 

n 3 = A *4 ( “n + k 2' Pn + l 2’*n + m 2- 6 n + n 2> ?n 

+ A) 

(F28) 


In the preceding fourth-order formulas given we have defined 

A “ ?n + l ' < F29 > 
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